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A General Inverse Problem for the Growth- Fragmentation Equation 
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Abstract 

Aggregation-fragmentation equations arise in many different contexts, ranging from cell division, 
protein polymerization, biopolymers, neurosciences etc. Direct observation of temporal dynamics 
being often difficult, it is of main interest to develop theoretical and numerical methods to recover 
reaction rates and parameters of the equation from indirect observation of the solution. Following 
the work done in [2] and [5] for the specific case of the cell division equation, we address here the 
general question of recovering the fragmentation rate of the equation from the observation of the 
time-asymptotic solution, when the fragmentation kernel and the growth rates arc fully general. 
We give both theoretical results and numerical methods, and discuss the remaining issues. 

keywords. Growth-Fragmentation equation ; general fragmentation kernels ; inverse problem ; 
eigenvalue problem. 
MSG Classification. 35Q92, 35R30, 45Q05. 

Introduction 

To model the behavior of a population where growth and division depend on a structuring quantity 
of the individuals such as size, the following mass-balance equation is currently used: 

( d d f+°° 

Trn{t, x) + c—{g{x)n{t, x)) + B{x)n{t, x) = 2 B{y)K{x, y)n{t, y) dy, t > 0, x > 0, 
Ct Ox Jo 

n{t = 0,x) = n°{x), x>0, 
_ g{x = 0)n{t, x = 0) =0,t>0, c> 0. 



Here, n denotes the density of the individuals structured by the size variable x at time t; the growth 
rate is given by g{x); the division rate B{y)n{x,y) represents the rate at which a given individual 
of size y gives birth to two individuals of size respectively x and y — x, whereas B{y) is the total 
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rate of division for individuals of size y. This physical interpretation of y) leads to the following 
assumptions: 

/•oo 

^(x, y) = Vx > y, / K(x,y)dx = l, K{x,y) = K{y - x,y). (2) 



By simple integration and symmetry, it leads to the following well-known relation, expressing the 
conservation of mass by the division process: 

oo y 

xk{x, y) dx = -. 

Problem ([T]) or its variants arises in many different contexts, ranging from cell division, protein poly- 
merization, telecommunication, neurosciences, and its mathematical study can provide useful infor- 
mation on the qualitative behavior of the phenomenon under consideration (see, among many others, 
[T2t [T3]). To be able to use it as a predictive model however, it is crucial to be able to estimate 
quantitatively its parameters (7, B and k. 

A first step consists in the use of the asymptotic behavior of this equation, as first proposed in 
[2]. Indeed, by general relative entropy principle it is proven (see e.g [111 [HI [5]) that under suitable 
assumptions on k, 5 and B one has 



00 
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with po = f nP{x)(j){x)dx and {Xo,N,(j)) is the unique eigenpair solution of the following problem: 



00 



+00 



c|-(5(x)iV(x)) + {B{x) + Ao)iV(x) = 2 / B{y)K{x, y)N{y) dy. 



gN{x = 0) = 0; N{x) > 0; / N{x) dx = 1, Aq > 0, 

cg{x) — ^{x) - {B{x) + XMx) = -2B{x) / ^(y, x)</.(y) dy. 
Ox Jo 

poo 

0(x)>O; / (/)(x)7V(x)dx = 1. 



(3) 



The use of this new problem allows to restrict the need for information to a non-temporal measure, and 
the problem becomes: How to recover information on y, B and k. from an experimental measurement 
of the asymptotic profile N and the global exponential rate of growth Aq of the population ? 

In the case when the equation models cell-division, direct measures of the growth rate y(x) is 
possible. Direct measures of k is also possible, by a study of the sizes of the two daughter cells born 
from a mother. The most delicate point is thus the measure of the division rate B, what implies to 
follow each cell from its birth to its division or death. 

In [2] and [3] , the problem of recovering the division rate from a measured was addressed in the 
case when the growth rate is constant, i.e g{x) = 1, and the daughter cells are twice smaller than 
their mother, i.e. when K{x,y) = 6x=y/2- In this case, Problem ([3]) writes: 

c-^N+{B{x) + Xo)N = 4.B{2x)N{2x). (4) 

Ox 



^Growth can naturally be balanced by death, by the addition for instance of a death term d{x)N{x) on the left-hand 
side of the equation. This would lead to possible nonpositive rates Aq, but our whole study would remain unchanged. 
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In this particular case, the inverse problem reads: How to recover H = BN{x) solution of 

C{H) = F{N), (5) 

with C : H ^ 4:H{2x) - H{x), and F{N) := cdj^N + AoiV ? The method used to solve Equation 
([5]) El strongly uses the analytical study of the operator C, and it was shown that the most efficient 
technique was then to view the problem as written in the variable y = 2x rather than in x (see the 
discussion in [3]). 

In this paper, we address the inverse problem of determining the cell division rate B when g and k 
are known - or guessed - functions, but fully general ; hence, we cannot apply anymore the inversion 
of the operator C as done in [3], and new tools have to be designed. 

We model the experimental measure of the distribution N by an approximation data Nf, of 
satisfying || — A^^^H < e for a suitable norm ||»|| H 

The paper is organized as follows. We first study the regularity of the direct problem, what is a 
necessary step for a better understanding of the inverse problem. In a second part, we investigate 
the inverse problem of determining B by the Quasi-reversibility and Filtering methods proposed in 
[2] and [3] and properly adapted to our general context. In a third part we develop new numerical 
approaches in order to recover the rate B following the two regularization methods ; we give some 
numerical illustrations of our methods. 

Main notations and assumptions 

We use the following notations. 

P := {/ > : 3/i,i^ > 0, lim supx ^f{x) < oo, liminf x'^/(x) > O}, (6) 

Ll : = {/, 3a > 0, / G LP(0,a)}, := L\^+,x^dx). (7) 

We work under the following technical assumptions, that guarantee well-posedness of Problem ([3]) 
as stated in [5] (we refer to that paper for a complete discussion and justification). 

y 

30<c<^, Vp>2, Dp:= j ^K{x,y)dx<c<^. (8) 



B G LlM*+) n P, 3 ao > 0, 5 G L^^(M+, a;°«dx) n P (9) 
y K compact in ]0, +oo[, 3 nik > : g{x) > \/x £ K (10) 



3 6 > 0, suppB = [b, +oo) (11) 

(12) 



3 C> 0, 7 > 0, 4^ G ; f k(z, y) dz < min (l, C(-)^) 

gix) Jo y ' 

B{x) . xB{x 



B(x) ^1 xB(x) , , 

e Ll ; —yi^ , ) +00. (13) 



^the method was first developed in [2], then investigated deeper and solved numerically in [3] in a deterministic 
setting, and in [7] in a statistical setting. It was also successfully applied to experimental data in 2]. 
''a more precise model for the measured data, in a statistical setting, can be found in [7]. 
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1 Regularity of the direct problem 



Before studying the inverse problem, it is necessary to have a proper knowledge of the direct one, 
which states as follows: What is the regularity of the map T : (c, B) — t- (Aq, N) solutions of Problem 
([3]) ? How can we define a proper definition domain for T ? 

In [3], Theorems 3.1. and 3.2 establish that the map Fq : B ^ {Xq,N) is Lipschitz-continuous for 
c = 1 fixed, g = 1, K = ^S^^y and division rates B such that < Bm < B < Bm < oo; in other 
words, for division rates uniformly positive and uniformly bounded. 

In this paper, we want to state such results for general growth rates g and division kernels k, with 
division rates B not necessarily uniformly bounded. Our study is thus first based on the well-posedness 
of this general eigenvalue problem ([3]), as performed in [5]. 

Let us first settle a proper definition space for the division rates B. Theorem 1 of [5] states that, 
under Assumptions ([2]) and ([8|)- (fT3|) . there exists a unique eigenpair (Aq, N, cp) solution of Problem ([3]). 
Hence, we first need that g and k satisfy Assumptions ([2]), ([8]l- (fT0]) . Then, to study the regularity of 
the map T : {c,B) — )■ {X,N), one needs not only that such division rates B satisfy Assumptions ([9]), 
(jlip and (|13p but also that they satisfy them uniformly. This leads to the following definition. 



Definition 1.1 Let g, k satisfying Assumptions ([2]), ([8]l- (fT0]l . For a constant b>0 and functions 
/o G in, /oo > oo, one defines the set 

^(^/o,/oo) := |i?GL;~(M;)nP, Supp{B) = [b<b,+^), j<fo, ^>/- 



In such a set, division rates B satisfy uniformly Assumption ()13p . what allows to use the powerful 
estimates proved in [5]. 

Under such assumptions, we also recall that we have the following results (see Theorem 1 in [5]) for 
the unique solution (Aq, A^, </>) to Problem ([3]): 

x'^gN G LP{R+) V a > -7, VI < p < +oo ; x'^gN G W^''^{R+) Va > 0, (14) 

and 

3k>0, C >0, (pix) <C{l + x''); 59^(/> G L°°(M+). (15) 

The two following fundamental estimates are straightfully obtained by integration on [0, 00 [ of Equation 
([3]) or ([3]) multiplied by x: 

r+00 

Ao = / B{x)N{x)dx, (16) 



JO 

X) g r+oo 

xN{x)dx = — g{x)N{x)dx. (17) 
/o -^0 Jo 

We are now ready to state our regularity result. 

Theorem 1.1 Let parameters g and k satisfy Assumptions ([8|)- ()10p . then 

i) The map T : {c,B) 1 — > {Xq,N) is: 

• continuous in {c,B) under the L°°— weak-* topology for B from any set x I?(6, /o,/oo) to 
X nL°°(M+). 
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injective. 



ii) Let moreover g satisfy ^ E Lq with 7 defined in Assumption ()12p . Then the map T is Lipschitz- 
continuous under the strong topology ofW^ x H /o, /oo)- More precisely, denoting 6B = 
B — B, A = \\B — S||l2(]k^), 6c = c — c, 6N = N — N, 6X = Xq — Aq, we have the following 
estimates, for k as in (I15|) ; 

\6X\<Ci{B,B)A + C2{B,B)\6c\, \\6N\\l2^^^) < C^{B,B)A, 

with 

^1 = ^"TT^II^^^Ti^^^' ^2 = ^^\\{l + x')gN\\LHR+) + m + x')gN\\L4B\\L^, 
where C > is an absolute constant. 
Proof. 



i) The continuity of the map T directly fohows from the proof given in [5], Theorem 1. Therefore, 
we only sketch the main steps and let the reader refer to this article. 

Let c„ — >• c > in and Bn B in L°°(M+). Denoting {Xn,Nn) the respective eigenpairs 
solutions of Problem ([3]) settled for {cn,Bn), we can prove the same uniform estimates for Nn 
as in [5] due to the fact that since Bn G 25(6, /o, /oo). Assumption is uniformly verified. 
Such estimates give strong compactness in for Nn, and hence, up to a subsequence, we have 
a strong convergence of {Xn,Nn) to (A, A^). Similarly, we prove A > 0, and passing to the limit 
in the equations for Nn, we deduce that (A, N) has to be the solution of Problem ([3]) settled for 
{c,B). Since such a solution is unique, the whole sequence {Xn,Nn) converges to it. 

Let us show by contradiction that F is an injection function. 

Let Bi G Lj^^(]R^) and q positive constants Vi € {1,2} such that {ci,Bi) ^ (c2,i?2) and 
T{ci,B,)=r{c2,B2) = {Xo,N). 

We then integrate the two equations satisfied by (Xq^N) against the weight x, to obtain 

' Aq Aq 

g{x)N{x) dx = — yN{y) dy = — yN{y) dy, 

what implies ci = C2- By the contradiction assumption we get Bi 7^ B2, so by making the 
difference between the following Equations (fTH]) . (fT9|) with consideration to the equality ci = C2 

d l"^ 

cr — {g{x)N{x)) + {Bi{x) + XQ)N{x) = 2 Bi{y)K{x,y)N{y) dy, (18) 

Ox Jo 

d 

C2^{g{x)N{x)) + {B2{x) + Xo)N{x) = 2 B2{y)K{x,y)N{y) dy, (19) 

Ox Jo 

we obtain after multiplying hy x^, p > 2 the following relation 

/•OD 

xP\Bi - B2\{x)N{x) < 2 / x^'l^i - B2\{y)N{y)K{x,y)dy. 

Jo 
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We integrate this relation on (0, oo) and due to Assumption ^ for p > 2 we deduce the fohowing 
strict inequahty : 

/•oo /"OO 

/ xP\Bi- B2\{x)Nix)dx < xP\Bi- B2\{x)N{x)dx, Vp > 2 
Jo Jo 

what is contradictory. 



ii) First, the fact that — G Lq impUes that for ah p > 0, G ((1 + x^) dx) . Indeed, by ([14]) and 

/oo f 
N'^{x){1 + x^)^ dx < oo, and it only remains to bound / A^(l + x^)^ dx. This 



smce 



X^T 

is given by writing A^^(x) = {N'^g^x '^^){ — product of an function with a Lq function. 
By making the sum between the two following equations 

(t){x)c—{g{x)N{x)) + (t){x){B{x) + Ao)iV(x) = 0(x)2 / B{y)K{x, y)N{y) dy 

Ox Jo 

N{x)cg{x)^(t){x) - N{x){B{x) + \Q)(t){x) = -2N{x)B{x) f x)(/>(y) dy 







we obtain 

d - d - - - - r°° - 

(5c</.— (5iV) + — (cgiV</.)(x) + ((/>iV[5 + Ao-i?-Ao])(x) = 20(x) / B{y)K{x,y)N{y)dy 

/•X 

-2N{x)B{x) / K{y,x)(t>{y)dy 
Jo 

we then integrate this equation on [0, oo) that leads 

6c (j)—{gN)dx + 6X cpN dx + (^N6Bdx = 2 5B{y)N{y)[ j (t){x)K{x,y) dx] dy. 
Ox Jo Jo Jo \Jo J 

So 

5\ j (/)(x)iV(x) dx = j 5B{x)N{x)y2 j (l){y)K{y,x) dy - (l){x)\ dx + 5c j gN—cpdx. 

The first term of the left-hand side gives the term with Ci{B, B) of the estimate on 5\ by using 
the fact that 3 A; > 0, ^ ^ ^ G -L°°(IR+). For the second term, we use the equation for (j) ^-iid 
write 

5c j gN^(pdx =^-^ j gN({B + Xo)(t> - 2B £ K{y, x)(t>{y) dy^ dx, 
and it provides the term with C2{B,B) in the estimate for 5X. 
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To prove the estimate on 5N, we make the difference between the two following equations 

c—{gix)N{x)) + {B{x) + Xo)N{x) = 2 / B{y)Kix, y)N{y) dy 
Ox Jo 

c—{g{x)N{x)) + {B{x) + Ao)A^(x) = 2 / B{y)K{x, y)N{y) dy 

Ox Jo 



5c—{gN)+c—{g5N) + [ {\o + B)N - {\o + B)N ± {\o + B)N 



we obtain 

dx d: 

That implies 



(^BN - BN±BN^{y)K{x,y)dy. 



Sc—{gN) + c—{g6N) + (Aq + B)6N 

Ox Ox 



2 / N{y)K{x, y)5B{y) dy - {5\ + 6B)N 
Jo 

l-OO 

+2/ B{y)5N{y)K{x,y)dy. 
Jo 



We recast the previous equation as follows 
d 

c—{g{x)6N{x)) + {Xo + B{x))6N{x) = 2 B{y)5N{y)K{x,y) dy + 5R{x), (20) 
Ox Jo 

with 

6R{x) = 2 / N{y)K{x, y)5B{y) dy - {6X + 6B)N - 6c— (gN) (21) 

JO Ox 

We can bound ||5i?(x)||j;^2 as we previously bound \6X\. The estimate on ||5A^||^2 thus follows 
from the following lemma. 



Lemma 1.1 Under the assumptions of Theorem M . l\ for a and k, with 5N defined as in Theorem \l.l\ 

and 6R defined by (j2ip . there exists v{c,B) > a constant depending only on the eigenvalue problem 
([3]) stated for given parameters c > and B £ L"^ (1 'D{b, fg, f^^) such that, for all c > cq > and 
B £ L'^ n V{b, fo, foo), one has 

i'll'^^llL2(M+) < II<^^IIl2(M+)- 

Proof. 

We argue by contradiction and assume that for a sequence > cq > 0, Bj^ £ L'^ Ci 'D{b, fo, foo), one 
has, for a vanishing sequence i^k, 

'^k\\SNk\\L2(R^) > \\SRk\\L2{R+), 

with 6Nk = Nk — N, Nk solution of Problem ([3]) stated for Ck and Bk, 6Rk defined by (f2T]) stated for 
6Nk. 

As for the proof of continuity above, compactness arguments as done in [5] lead us to extract a 
converging subsequence Nk — ^ N strongly in L^, so 5Nk — )• 5N strongly in L^. Moreover, estimates as 

in [5] imply that Nk is uniformly bounded in (we write iV| = x^'^^g'^N'^ — ^ and use the assumption 
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— € together with the result (I14p . result which is uniform for all Nk), hence 6Nk satisfy Equation 

g 

()20p with ||(5i?fc||^2 — )• 0. Passing to the limit, it implies that 5N satisfies Equation so by uniqueness 
of a solution we have 6N = CN for a given constant C G M. Since J Ndx = J N^dx = 1, we have 
J SNdx = C = : it is contradictory with our assumption on (i^k)- □ 



2 The inverse problem and its regularization 

As in [U [7] , we consider the problem of recovering the cell division rate B and the constant c from 
the a priori knowledge of the shape of the growth rate g{x) and the experimental measure of the 
asymptotic distribution and exponential growth Aq- To model this, we suppose that we have two 
given measurements Ap G L^nL°°(M+) and > such that ||A — A^IL,/ n , \ < lAn — A^l < eu 

& fc \ - II '''' L^(^{l+xP)dxj — ' I u fci _ i_i 

The problem is: How to get estimates {c£,B^) of {c,B) solutions of 

c—{g{x)N{x)) + {B{x) + \o)N{x) = 2 / k(x, y)B{y)N{y) dy. (22) 

Ox Jo 

First, it is clear that B cannot be recovered from Equation (j22p when the distribution A vanishes: 
our inverse problem consists in recovering H = BN rather than B directly. Our problem can now be 
viewed as: How to recover (c, H) solution of 

C^{H){x) := H{x) - 2 / k(x, y)H{y) dy = -c—{g{x)N{x)) - AoA(x) (23) 

Jo Ox 

when we have measurements (A^, A^) of (Aq, A^)? 

Secondly, since the measure Ag is supposed to be in L^, there is no way of directly controlling 

d 

— {gN^) even if g is known (see Section 2 of [2] for a discussion, or yet [l]).To overcome this difficulty. 

Ox 

two regularization methods were proposed in [21 [3] for the particular case of division into two equal 
cells, i.e. when K{x,y) = 5^_j^/2) ^ third method has also been proposed in [11], and a statistical 
treatment to estimate the derivative in [3- Indeed, looking at the problem in terms of H = BN and 
not in terms of B makes it almost linear in H; almost, because Aq being also measured, the term AqA^ 
can be viewed as quadratic. Hence, the classical tools designed to regularize linear inverse problems 
(see [U) can be used, as illustrated by the three foreseen methods, as soon as the operator can be 
inverted. 

This is the third and last difficulty: inverse the operator defined by Equation (j23p . None of the 
three regularization methods of [21 [3l [H] can be directly applied here: indeed, they strongly used the 
fact that for the kernel k = S^^r, the left-hand side of Equation (j23p simplifies in 4BN(2x) — B{x), 
and can be viewed as an equation written iny = 2x. Then, a central point of the proofs in [2] as well as 
in |3] or [11] is the use of the Lax-Milgram theorem for the coercitive operator C : H ^ ^H{y) — H{^). 

Nothing such as that can be written here, and the main difficulty, numerically as well as theoretically, 
is to deal with a nonlocal kernel / K{x,y)H{y)dy. The operator C is replaced by For k = S^^r, 
Ck has been proved in [3] (Proposition A.l. in the appendix) to be coercitive in L?'{x^dx) if p > 3, or 
in contrary C is coercitive if p < 3. Due to the nonlocal character of the kernel, it seems more natural 

*See [7] for a statistical viewpoint on the data (A'e,Ae) : supposing that A^g £ means that we deal with some 
preprocessed data. However, once the problem is solved in a deterministic setting, as we do in this article, it is immediate 
to apply the method of [7| to this general case. 
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now to look for cases when the first part of the operator C^-, i-e. identity, dominates the nonlocal part 

oo 

2 J u{y) n{x , y) dy . This is expressed by the following proposition. 

X 

Proposition 2.1 Let k satisfy Assumption ([2]). For r,q > 0, we define the following quantities: 

oo y 
:= sup I —K{x,y)dy, := sup —K{x,y)dx. (24) 

X J y y J y 

X 

If < r < p are such that 

CrDp_r < I, (25) 

Then for all F G L^(M+,x^dx) there exists a unique solution u G L^(x^d2;) to the following problem: 

oo 

u{x)-2 J u{y)K{x,y)dy = F, (26) 

X 

and we have the following estimate 

IWWl^ixPdx) < =11-^1 lL2(xPdx)- 

Proof. We define the bilinear form 

r+oo r+oo r+oo 

A{u,v) = / u{x)v{x)x^ dx — 2 / u{x) / K{x,y)v{y) dy x^ dx =< u,v >x,2(^xp dx) 

Jo Jo Jx 

where <, >^2 denotes the scalar product. We apply the Lax-Milgram theorem in L^(x^dx). Indeed, 
we have by Cauchy-Schwartz, for any constant C > : 

roo roo roc /•y 

B{u,v) = 2 / xPu{x)v{y)K{x,y) dy dx = 2 / / x^u{x)v{y)K{x,y) dx dy 
Jo Jx Jo Jo 

oo y 

f f x^ 1 xP~^ 

< J j {Cu^{x)x'P— + —u^{y)yP--^)K{x,y)dxdy 



oo cxD oo 

/r rj^r If f x^~^ 

u^{x)xP{ J —K{x,y)dy)dx + — J u'^{y)yP{ J -^K{x,y)dx)dy 



oo 

< {CCr + %^) y" u^{x)x^dx, 





The minimum is reached for C = \/ jj^^ So finally we have 



what proves the continuity of the biUnear forms B and A. Moreover, it imphes 

A{U,U) > (1 - 2y/CrDp_r)\\u\\l2^^r>dx) = ^Ml^xP dx^ 



with /3 = 1 — 2y/CrDp^r > under assumption (j25|) . It ends the proof of Proposition 12. 1[ 

Remark 1 Assumption ()25p can be linked to Assumption ([8|). For self-similar kernels n{x,y) = 

-A;o(— ) with / kQ{z)dz = 1, defining Ir = J z^'kQ{z)dz we obtain Cr = Ir~i <ind Dq = Ig, so that 
y y Jo 

Assumption ()25p is reduced to 

If— lip— r ^ ~- 

For the equal mitosis kernel K,{x,y) = 6^^y, since Iq = 2'^, Assumption ()25p is verified for p > 3 : we 
recover part of the result of the proposition of f^. It corresponds to the cases when the first part of 

the bilinear form (i.e., J uvx^ dx) dominates the second one ( J J K(x,y)u{x)v{y)x^ dx dy). 

For the uniform kernel kQ{z) = 1, the equality Iq = here again implies that Assumption (j25p is 
verified for p > 3. 

More generally, for these self-similar kernels Assumption ([8]) is equivalent to Ir>2 < ^- Since Iq = 1 
and Ii = ^ by Assumption taking r = 2 in Assumption ([25]) leads to 

and due to Assumption ([5]) this is always true for p > 4. 
2.1 Filtering method 

This regularization method consists in looking for a solution i^e.a of the following regularized problem 

r+oo ( d \ 

C^{He,o){x) := He,a{x) - 2 y)He,aiy) dy = * i-ce,aj- {gix)N,{x)) - XsN^ix) j , (27) 

where Pa is a mollifiers sequence defined by 

1 X f'°° 
Paix) = -pi-), p G C~(M), / p(x) dx = 1, p > 0, Suppip) C [0, 1]. (28) 
a a Jq 

One notices that c^^q is uniquely defined: indeed, integrating Equation ([27|) against the weight x leads 
to 



xN^ dx 

• (29) 

Pa * {gNe)dx 



We want to study the well-posedness of this problem and estimate the distance between Bir a = jf^ 
and B in order to choose an optimal approximation rate a. This is given by the following result. 
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Theorem 2.1 Let g, B and k satisfy Assumptions ^ and ([5|)- (|13|) . and moreover — ^ Lq with 7 

defined in Assumption (jl2p . Let (Aq, N) the unique eigenpair solution of Problem ^ (as stated in f^). 
Letp>l satisfy Assumption ([25]). Let Ne e L'^ n L°°{R+) and Xg > satisfy \\g{N - Ne)\\L2^^p ^^-j < 
e\\qN\\T2(„p A^\, lAe — Anl < eAn, II-/V — -/VglL,/, , \ <e||A^|L-,/', , \ and 

liy ll_L (x'^nxjj \ e ui _ U) II '^''L^\^{l+x+g{x))dx) — " [(l+x+g(x)) dx) 

Then there exists a unique solution H^^a £ L'^{x^dx) to Problem (|27p . 
Defining Bg^^a '■= Xn^ aix)^0'^e,a/Ne,a we have the following estimates: 

\ce,a-c\<C{p,p,N){a + e), (30) 
\\Bs,a - B\\L2(^.-,pN2dx) < C{p,p,N){a + (31) 

where C is a constant depending onp, moments of p and Aq, lb^llj:^2 ((i+^p) d^) ' ^^^h^{(i+x+g(x))dx) ' 
|biV||H.M(dx) and l|A^II^,((,.,+i)dx)- 

The estimate (I35p of Theorem 12.11 relies, on the one hand, on the estimate of Proposition 12.11 and, 
on the other hand, on general approximation properties of the mollifiers, as expressed by Lemma |2. II 
right above. 

Lemma 2.1 Let p > 1, / G L^((x^ + 1) dx), pa a mollifiers sequence defined by (|28l) andO < a < I. 
Then we have the following estimates. 

1- 11/ * PaWl'^lxP dx) 

<C{p,p)\\f\\ 

L^(^[xP+i)dx)' ^iP^P) "^^y depending onp and moments of p, 

d 1 

2. \\-Q^{f*Pa)\\L^{xPdx) < ~^(P^P)\\f\\L2(^{xP+i)dx)' ^^if^C{p,p) only depending on p and moments 
of p and ■§^p. 

3. ||/*Pa-/||L2(.Pd.) <C(p)a||/||^i(,Pdx) 2//Gi/H(l+^^)d^) 
^. ll/*Pa-/llLi <C(p)a||/||H.M 

5. ||Pa*/||Li<ll/llLi- 

Proof. The proof of this result is classical and relies on Minkowski inequality for convolution products 
; we let it to the reader. 
Proof of Theorem 12.11 We decompose the left-hand side of Estimate (j35p as follows 

\\^e,Ci^ - B^^L^ixP dx) = ||^£,a(^ — + - -^£,a) + -f^e.a - -BAf||j^2(a,pdi) 

< ||-S£,a||L°° ^11^^ - ^(Al?-{^xPdx) + ll-^a " ^ e,a\\l?'{xP d^ + ll-^^e," " ^^^L^ixP 



On the right-hand side, the first term is bounded by C(p, /3)a|| A'"!! ^, /, , , \ due to Lemma [2.11 
Estimate O The second term is bounded by C(p, p)e||A^|| / \ due to Lemma [2.H Estimate 

[U applied to f = N — N^. For the third term, we apply Proposition 12.11 to u = He,a — BN and 
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F = Pa * (^e,a-^ {gNe^ + XeN^ — ^c— (^g{x)N{x)) + AqA^^ . We treat these terms in a similar 
manner. Let us detail briefly the most binding term (with the notation Lp = L'^{x^ dx)): 

\\pa* Ce,a-g-{9Ne) - c—{gN) \\l2 < c,,«M|— * (gN^ -gN)\\L2 + \\pa* g-igN) - ^{gN) ^2] 

d 

+ \Ce,a - c\\\-^^{gN)\\L2. 

The first term is bounded by C—\\gN\\i^2(^r^p^x-^ by Lemma lXTl Estimatel^ the second one by Ca\\gN\\^2(^xp 
by Estimate [3l For the third term we write 

J xNfrdx ^ J xN{x)dx 
f Pa* gNe dx / g{x)N{x) dx 
I ^ / xNi; ^ ^ J xN dx ^ ^ J xN dx ^ f xN dx 

J Pa* gNe dx " j Pa* g^e dx ^ J Pa * g^e dx ^ J gN dx 

f xN^ J xN dx J xN dx J xN dx 

< \Xe— — — ; — — I + lAe-;; — — ; Aq" 



J Pa* gNe dx f Pa* gNe dx f Pa * gNe dx f Pa * gNe dx 

J xN dx J xN dx 

'jPa*gNedx " ^°iVVd^'- 

The assumptions of Theorem 12.11 together with Estimates H] and [5] of Lemma 12.11 give the estimate for 
\ce,a — c| and ends the proof. 

2.2 Quasi-Reversibility Method 

To regularize the exact inverse problem (j23p . the so called quasi-reversibility method proposed in [2] 
for the case k = dx=^ consisted in adding a term derivative a-^{BN(2x)) with a small a > to 
the right-hand side of Equation ()22p . viewed as an equation taken in the variable y = 2x. The main 
difference is that we need here to take this term in the variable x and not 2x due to the general form 
of the nonlocal kernel k. We choose to define, for a > and A: G M, the following regularised problem 

Cl{He){x) := ax-^ — {x''+^He{x)) + He{x) - 2 J K{x,y)He{y)dy = -c,,, — (giV,(x)) -A,A^,(rE), 

He{0) = 0; 0<a<l, p£R. 

(32) 

This equation has to be understood in a distribution sense in M-|_ undowed with the measure x^ dx. We 
moreover assume that Supp{Ne) C K^. Other adaptations would be possible, all consisting in adding 
a small term derivative of the form zizafi{x)-^{f2{x)BN{x)), with a > and a boundary condition 
taken either in x = if a > or x = +00 if q < 0. Numerically indeed, a < proved to give better 
results (see below Section [3?3l) . The key point is to check that the regularized operator £^ satisfies 
Proposition 12.21 below. 

The choice of Ca,e is not directly given by integration of the equation, contrarily to the case of [2]. 
Neglecting the regularisation terms involving a, we thus define, as for the exact equation ([22]): 



Xe J xNe{x) dx 

Ca.e = f • ("^3) 

g{x)Ne{x)dx 
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Theorem 2.2 Let g, B and k satisfy Assumptions ([2]) and ([8])- (jl3p . and moreover ^ € Lq with 7 



g 

defined in Assumption ()12p . Let {\o,N) the unique eigenpair solution of Problem ([3]) (as stated in fBjj). 
Letp>2 satisfy Assumption (gSD- Let £ L^nL°°(]R+), Supp{Ne) C M^, and > satisfy \\g{N- 

A^e)||L2(.Pdx) <e||5iV||L2(xPdx), |A,-Ao| <6Ao, ^^II^IIli d.) 

11^^ - ^llL2((.p+i)d.) ^ ^ll^llL2({xP+i)dx)- ^ dx) 6e so/uiion to Problem 

Defining -Be, a := Xn^{x)^qH^/N£ we have the following estimates: 

\ce,a-c\<C{p,N)e, (34) 

\\Be,a - -B||L2(xPAr|dx) < C{p,N){a + (35) 
where C is a constant depending on p, k, Aq, \\'^^\\fji^(^i_^_^p+i-^^^y ll"^llLi((i+x+g(a:))dx) ' 

\\gN\\Li and 

\\N\\ I \. 

" "/fi ((xP+l)dx] 

Proof. The estimate for \c^^a — c| is obtained in a similar manner as for the filtering method. 
For the estimate for we first write 

||-Be,QiVe - BN^^iji^^vAx) < " Li^\\L'^{xv Ax) + 11^^ " IIl2(xp dx) • 

The second term of the right-hand side is simply bounded by 112,2(2.^^3,) \ N e\i?'(xT Ax) ^ £||-S||L2(3;Pdx)- 

For the first term of the right-hand side, as for the filtering method, we decompose Li^ — BN, 
and for this we need to establish some regularity properties of the operator defined in Equation 
(j32p and designed to approximate £k- This is given by the following proposition, which is for the 
quasi-reversibility method the equivalent of Lemma l2.1l for the filtering method. 

Proposition 2.2 Letp>2, F = fi + ^/a with fi G L^ {{1 + x) dx) n L"^ (xP dx) and f2 G H^{{1 + 
xP)dx) n W^'^{x dx). Let K,g,p satisfy the assumptions of Theorem \2.2l There exists u G L^{xdx) 
solution of the following problem, where k ^ —2 and < a < 1 : 

C'^{u){x) := ax~''-^{x''+\) + C^{u) = F. (36) 
Ox 

Moreover, we have the following estimates for a constant C > only depending on g, k, k and p : 

1- \W\\l'2{xPAx) < r-^ WL'WL^ixPAx)^ 

1 — Z a/ Op 

C 1 

^- IkllLaCxPda;) < — ll/l + /2(1 + -^)\\l^(xP Ax)- 

Proof. Let us first establish the existence of a solution in L^{xdx). We rewrite (|36p as follows 

ax—{u{x)) + {a{k + 1) + 1)^(3;) = 2 j k{x, y)u{y) dy + F{x). ^^^^ 
u(0) = 0, p > 2. 

We consider v G L^(]R+, xdx) and define u = T{v) the explicit solution of 

d 

ax—{u{x)) + {a{k + 1) + l)u{x) = 2 1 v{y)K{x, y) dy - F{x), 

Ox Jo 
n(0) =0, p > 2. 
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Let vi and V2 two functions of L^{W^,xdx) associated to ui and U2 then by doing the difference 
between the two equations satisfied in the one hand by ui , vi and in the other hand by U2 , V2 we have 

d f°° 
ax^{6u{x)) + {a{k+l) + l)6u{x) = 2 / 6v{y)K{x,y) dy, with. 6u = U1 — U2 and 6v = V1 — V2, (38) 



dx 

what implies the inequality (see |13j . prop. 6. 3 for instance) 

d f°° 
ax^\5u{x)\ + {a{k + 1) + l)\5u{x)\ < 2 / \5v{y)\K{x, y) dy. 

Cx Jo 

Multiplying by x and integrating on [0, 00 [ we deduce the estimate 

x\6u{x)\dx < — . / y\5v{y)\dy. 

a{k - 1) + 1 Jo 

This proves that T is a Lipschitz function and we deduce the existence of a solution u £ L^(M_|_,xdx) 
by the Schauder fixed point theorem. 

For the first estimate, we multiply Equation (I36p by x^u and integrate from to x. Using that 

o ft 2 

uxP-^-^{x''+\) = {k + l)xV + xP+^-^i'^), it gives 

rx i-x rx 

a{k + l)xPu^{x)dx + -xP'^^u^{x) + I Ci^{u){x)u{x)x'p dx = I F{x)u{x)xP dx. 
2 Jo Jo 



From this, we deduce 



Ck{u){x)u{x)xP dx < J F{x)u{x)xP dx. (39) 

Applying the coercitivity on L^(.^ dx) of the bihnear form ^(n, .) = / £.(.)..^ dx we get immedi- 
ately the first estimate. 



For the second one, we integrate by part, on the right-hand side of Equation (j36p . the term with 
^/2, and use the equation to express -^{u) with the other terms of the equation: 

f d f d f p f d 

j {-Q- f2)uxP dx = -J f2-^{xPu)dx = - J -f2UxPdx- J f2xP—{u)dx 

= - [ ^f2UxPdx+ [ 'l±lf2UxPdx+ [ J-f^{C^{u){x)-h{x)--^f2)xPdx 

J X J X J ax ^ Ox 

= J ^^l~ ^ f2UxPdx + ^(^1 ^/2(£.(n)(x) - Mx))xP dx + ^ I flxP-^dx^ 
C p 

< —\\f2\\L^{{xP~'^+xP)dx){\\u\\L^(xPdx) + ll/l||L2(xPdx)) + ~\\h\\ l2(xP-2 dx) 

Together with the first estimate, it provides the desired inequality. 
We are now ready for the proof of Theorem 12.21 We see that can be viewed as solution of 

d 

Equation (j36p with F = —Ca,e-^{g^e{x)) —^e^e{x), whereas H = BN would be solution of ([36]) 

Ox 
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if a = and F = —c-g-[gN{x))—XoN{x). To isolate in the error term the contribution due to the 

a— regularization from the one due to the measurement error e, we thus define an intermediate function 

d 

Ha as the solution of Equation with F = —c—(^gN{x))—XoN{x). We then write: 

Ox 

The function — Ha is solution of Equation ()36p with 

d 

Fe,a = -Q-{- Ca,egNe{x) + CgN (x)) - XeN,{x) + AoA^(x), 



X 



and we can use Estimate [2] of Proposition 12.21 to obtain an error term in the order of — . The difference 

a 

d 

Ha - H is solution of Equation ([36]) with F = -ax-^ — {x^^'^BN), and we can use Estimate [J of 

Ox 

Proposition 12.21 to obtain an error term bounded by ^^\\BN\\^-^^^_^^^^^^^^^. It ends the proof of 
Theorem 12. 2i 



3 Numerical approach of the inverse problem 
3.1 The direct problem 

Assuming that the division rate B, the growth rate g and c > are known, we solve the time-dependent 
problem ([1]) and look for a steady dynamics. 

We choose to split the time evolution of the problem into its conservative advection part and into its 
gain and lost part by division as follows 

d d 

—n{t, x) + c—{g{x)7i{t, x)) = 
Ot Ox 

d [°° 

-^n{t, x) + B{x)n{t, x) = 2 B{y)n{x, y)n{t, y) dy. 
- Ot Jx 

We use an upwind finite volume method with computation length domain L and grid number points 
ka: Xi = iAx, < i < ka with Ax = L/ka 

nf = / n(k At, y)dv, / n(A;At + s, x,) ds « nf+^ 



1 f^i+ ^ 1 r 

— / ^ n{kAt,y)dy, — / n{kAt + s,Xi)d 
Jx, 1 At Jq 



For the time discretization one can choose, thanks to the CFL (Courant-Friedrichs-Lewy) stability 

condition, the time step At < g with the notation gi = g{iAx) and Bi = B{iAx). 

max {Bi + —gi) 

iel,...,ka Ax 

The numerical scheme is given for i = 1, ka hy Uq and 
• First for the conservative equation 

\k 



the interface fluxes (5'f^)fj_i/2 defined by upwind method. 
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Second for the gain and loss part by cellular division we compute 



n. 



1- AtBi)n'^^''' + 2AtJ^I: 



r^ka 

where J"/^ « / B{y)n^^^''^{y)K{xi,y) dy. 
• At last we renormalize the discrete solution by 



ka 



what allows to have h^^^ > N, > NiAx =1, iVj > 0, where is the dominant 

fc— >-co ^ — ' 
i=l 

eigenvector for the discrete problem associated to the following steady equation 
c—{g{x)N{x)) + {B{x) + \o)N{x) = 2 / B{y)K{x, y)N{y) dy 

Ox Jo 

with Aq the dominant eigenvalue associated to N. 

3.2 The inverse problem without regularization 

As illustrated in [3], solving numerically Equation (j23p without regularization is unstable. Indeed, 
this recovering naive method gives bad reconstructions of H = BN as soon as the observed A^ is 
irregular (see above the estimates and see also [3]). Here, what we moreover observe is that, at the 
neighborhood of x = 0, the solution explodes. As an example of this we consider a length domain 
L = 25 and the total number grid points ka = 300. By an upwind method we compute numerically 
the Equation (|23p and compare the result with that obtained by solving the direct problem (D.P). 

3.3 The inverse problem : Quasi- Reversibility discretization 

In this section we numerically investigate the regularization of the inverse problem ()23p by the Quasi- 
reversibility method. It is based on Equation (j32p . that we rewrite, dropping the index e, as follows 

-ax-^ — {x''+^B^[x)N{x)) + c—{g{x)N{x)) + {B^{x) + Xo)N{x) = 2 / k{x, y)B^{y)N{y) dy. 

Ox Ox J X 

(Bc.A)(0)=0, {B^N){oo) = Q- < a; /c G 1R+. 

Assuming that N and A are measured, we first define c by (|33p and then look for an estimate of the 
division rate B^. For this, we put the notation 

d 

= Bo,N and L = -c—(gN) - XqN. 

Ox 
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Figure 1: Numerical reconstruction of BN by "brute force method" with the choice g{x) = x and 

l^{x,y) = -lx<y 

y 



By a standard upwind method we obtain, when dropping the index a, the following discretization 

^k + l TT + 1 TT \ 



j=i 



Ax 



with L 

Ho = and Hi = 0, V / > ka 
By developing this discrete equation we obtain 



\QNi-c{ — ), \li = l,...,ka 



-a- 



ka 



2Kj_j+iAx j ffj+i + I 1 + m — 2ft;j^jAx j i^i — 2 HjKijAx = Li, V i = 1, fea. 
^ ^ ^ j=i+2 

We rewrite it under matrix shape A x H = L with A a matrix of coefficients of size ka x ka; H is the 
unknown vector of size ka and L is a known vector of size ka. 

The matrix A being a upper triangular one, we can solve directly the linear system thanks to the 
following iterations 



H, 



ka 



A 



ka,ka 



Hi = -^ihi - ^ Ai^jHA; y i = ka-l,...,l 



j=i+l 

The matrix A satisfying Aij > for j > i + 1, we can choose Ax small enough so that Ai. 
1 + ai — 2k j jAx > for all i. This guarantees that no oscillations appear. 
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3.4 The inverse problem: Filtering discretization 



This section is devoted to the numerical discretization of the inverse problem (j23p by the Filtering 
method based on Equations ()27p -(f28]l. The aim is to numerically solve the Equation (I27p that we 
rewrite when dropping the index e as follows 



xPcag-{g{x)NUx)) + xP{B^{x) + Xa)N^{x) = 2xP K{x,y)B^{y)N^{y) dy, 



(5«iV«)(0) = 0, a>0, 
with Na = N * Pa and pa a sequence of moUifiers. 

As previously, we want to estimate from a measured density N and Malthus parameter A. We first 
define c by (j29p . We then rewrite the regularised equation as follows 

roc Q 

xPBa(,x)Na{x) - 2 / xPK{x,y)B{y)Na{y)dy = -xfc, — (5(x)iV,(x)) - X^xPN^ix) 

For the convolution terms arising in the previous equation we use the combination of the Fast Fourier 
Transform and its inverse which we respectively note by F and F* then we define the mollifiers pa by 

its Fourier transform: PniC) = — , 

This leads to the following approximations 

Na ~ F*(p«(OF(iV)(0); §-j9Na) ~ clGN^ = F* {i^pa{0F{9N){0) . 
For the discretization we put the notation 

F[a = BaNa and La = —CadGNa — Xa^a 



then in each grid point Xi = iAx we obtain when dropping the index a: 
( Ho = 

ka 

x^{l — 2Ki^iAx)Hi — 2 x^HjKijAx = x^Li ; V i = 1, ka. 
j=i+i 



We rewrite this previous discrete equation under matrix shape A x H = L with A the matrix of 
coeficients which is an upper triangular one and of size ka x ka. 

The shape of the matrix A allows to use adequately the LU iterative numerical method, and then we 
deduce the following iteration 



L 



ka 



ka 



Aka. 



ka 



i = ka — 1, 1 



j=i+i 
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4 Numerical Tests 

For the numerical tests we use as input data the noisy one which correspond for e = to the 
eigenfunction N obtained by solving numerically the long time behavior of the direct problem in 
section 13. li The direct problem is solved in the length domain L = 25 for ka = 300 number grid 
points with two differents initial data: a step initial data and a maxwellian one, as follows 



and the steady solution is taken when ||n(t,x) — N{x)\\ii <= 10 

In order to show the unique asymptotic profile of the direct problem we plot in pictures Figl2] the 
steady cellular density N related to the two previous initial data with different values of c, B and with 




(40) 



the choice g{x) 



= X 



2 and K{x,y) = -I{x<y}- 
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X X 

Figure 2: Direct problem g = x^^^: Top left: Step initial function. Top right :Maxwellian initial 
function. Down left: Steady solutions of cellular density with c = 1. Down right: Steady solutions of 
cellular density with c = 0.5 . 



4.1 Numerical reconstruction of BN in the noiseless case £ = 

For the case where the input data are exactly known i.e. for e = 0, we recover thanks to the Quasi- 
reversibility and Filtering methods the division rate B by computing numerically the value of BN with 
N obtained by solving the direct problem with high precision and for various choices of the division 
rate B as shows in figure Figl3] below. 
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Figure 3: Various choices of B to solve the direct problem 



In figure FigJSl the Tray-function is defined as follow 

0, for X <2 



B{x) = < 



^ , for X G [2,15] 
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1, for X > 15. 

Then with the notation D.P for the direct problem we obtain: 

• For the choice k(x, y) = — Iix<:y\ 

V 



25 
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Figure 4: Numerical reconstruction of BN for each regularization method in the case e = 0. 
Top left :g{x) = x. Top right :g{x) = x^/^. Down ■.g{x) = x^l"^ . 
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We measure the relative error in norm by 



error 



\BN - {BN)e,4^, 



\BN\ 



(41) 



L2 



where BN is the exact numerical solution of the direct problem and {BN)^^a represents the numerical 
reconstruction either by the Quasi-reversibility method or by the Filtering one. So we obtain for 
instance for the given parameters g{x) = x, c = 0.015 and k = 2.34 the following reconstruction error 
of the division rate function of a. 
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Figure 5: Numerical errors for e = with different choices of B in the direct problem. Left : errors 
by Quasi-reversibility method. Right : errors by Filtering method. 



For the choice k{x, y) 
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Figure 6: Numerical reconstruction of BN for each regularization method in the case e = 0. 
Top left:g{x) = x. Top right: g{x) = x^^^. Down : g = x^/^ . 
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We measure the reconstruction error thanks to the relation (j4ip for the given parameters 
g{x) = x^/^, c = 0.5, k = 2.34 and we obtain the following representations 
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Figure 7: Numerical errors for e = with different choices of B and c in the direct problem. 
Left : errors by Quasi-reversibility method {k = 2.34) . Right : errors by Filtering method. 



4.2 Numerical reconstruction of BN in the noisy case e ^ 

For this case, we consider as input data the values of the solution of the direct problem in which 
we add a multiplicative random noise uniformely distributed in |] (see [7] for a more precise 
statistical setting of noisy informations). The nonnegativity of the data is insured by the choice 



Ne = max{N{l + le),0), h^,^], e£ [0,1]. 



Then with these noisy data we numerically obtain 



For the case k{x, y) = —I^x<y} 
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Figure 8: Numerical reconstruction of BN by the measured data for different values of e with the 
choice B{x) = exp(-0.08(x - 12)^), c = 0.015 and g{x) = x. 
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For various choice of the parameter e we compute the relative error thanks to the relation (j4ip 
and obtain the following representations 




Figure 9: Numerical errors for different values of e / with B{x) = exp(-0.08(j; - 12)2), ^ ^ q g^^g 
and g{x) = x m. the direct problem. Left : errors by Quasi-reversibility method {k = 2.34) . Right : 
errors by Filtering method. 



Remark 2 Let us note that for data with high noise values i.e. e > 0.075 the regularization by 
Quasi-reversibility method gives numerically better results than the Filtering one which creates big 
oscillations. 

1 X 

• For the choice K(x,y) = -ko(—) with kq ~ AA(i i) 

y y 

Discussion 

As shown by the numerical illustrations above, and after that we tried many different shapes of 
regularization (trying for instance a wide variety of k and p, with zta, in the quasi-reversibility method) , 
our simulations still present some delicate points. Indeed, even if the regularization methods prove to 
give better result than the naive "brute force" method as shown by FigureEl the gain remains relatively 
small, and the regularizing parameter a has also to remain small to avoid wrong reconstructions. Due 
to this small regularization, as shown by Figures [TUl the noise is filtered but not as much as we 
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Figure 10: Numerical reconstruction of BN by the measured data for different values of e 7^ 
with the choice B = min(l, j^), c = 1 and g{x) = x^^^. 
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Figure 11: Numerical errors for different values of e 7^ with B = min(l, j^), c = 1 and g{x) = x^/^ in 
the direct problem. Left : errors by Quasi-reversibility method {k = 2.34) . Right : errors by Filtering 
method. 

hoped first - especially for smaller x, that are farer from the departing point of the algorithm. Finally, 
the parameter a needs to stay in a confidence interval, selected, for a given growth rate g{x), from a 
range of simulations carried out for various plausible birth rates (see for instance Figures [U [Tj [9]) . 

5 Conclusion 

We have addressed here the problem of recovering a birth rate S of a size-structured population from 
measurements of the time-asymptotic profile of its density, in the general case when a given individual 
can give birth to two daughters of inequal sizes. Compared to the work carried out in [21 [3l [11] this 
last assumption has raised new difficulties, the principal one being that we have no other choice than 
considering the equation from the "viewpoint" of the daughter cell - what implies to take into account 
the nonlocal integral term. We established theoretical estimates and built numerical methods to solve 
it. As shown above by our numerical illustrations however, some issues still remain to be solved, 
especially the behavior of the algorithm for smaller x and the cancellations of oscillations (also present 

in mm)- 
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